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We introduce an exact Monte Carlo approach to the statistics of discrete quantum systems which 
does not rely on the standard fragmentation of the imaginary time, or any small parameter. The 
method deals with discrete objects, kinks, representing virtual transitions at different moments of 
time. The global statistics of kinks is reproduced by explicit local procedures, the key one being 
based on the exact solution for the biased two-level system. 
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1. Quantum Monte Carlo (MC) simulation is the most powerful, if not the only available method of getting exact 
results for complex systems Jl|-0] , where analytic solutions are not possible and exact diagonalization methods do not 
work because of the enormous Hilbert space. One conventionally starts from the general expression for the thermal 
average (A) = TrAe~P H /Z(/3), where Z(/3) = Tre _/3ff is the partition function for the Hamiltonian H, and converts 
the trace into the Feynman path integral gj. To sum over all trajectories with the weighting function e~ s defined 
^**^ , by the action S the imaginary time interval [0, 0\ is divided into Np = (3/ At smaller intervals of width At, and the 
actual trajectory is described by specifying the system state \at) at each time = kAr where k = 0, 1, Np. To 
ensure summation over trajectories with the near-optimal action, the complete sum is replaced by the stochastic sum 
(the known Metropolis algorithm jj))); in the most commonly used schemes the probabilities W\ t % to have the two 
trajectories {afc}i,2 in the sum are related as 
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ON ' W 1 /W 2 = e S2 ~ Sl . (1) 
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After some trajectory is accepted to the statistics, the MC process generates small variation to this trajectory which 
is cither accepted or rejected according to Eq. ([!]) . 

An approximate treatment of noncommuting operators originating from the artificial discretization of the time 
variable at points where the trajectory changes, i.e., when \otk) ^ |afc+i) ( in what follows we call such points 
"kinks"), results in an intrinsic error ~ At 2 (see, e.g., [f[o|). Taking smaller and smaller values of At eventually 
eliminates this intrinsic error, but then the algorithm starts accumulating statistics too slowly. Consider, as a typical 
example, the Hamiltonian of a particle on a lattice 



h = -tj2 A a J + Yl Uini ' ( 2 ) 
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where a\ creates a particle on site i, t is the hopping matrix element, rii = aja^, and < ij > conventionally denotes 
nearest neighbor sites. In this case the set of states {a} may be chosen to describe particle coordinate. The typical 
separation in time between two adjacent kinks is of order l/t and independent of At, so that for small At there 
are some l/(tAr) ^> 1 time intervals between the kinks. Now, if the MC procedure will suggest to create a new 
kink-antikink pair then the corresponding variation of the trajectory will be most probably rejected since the ratio 
W new /W id ~ (tAr) 2 is proportional to the square of a small parameter. On another hand, when the MC procedure 
suggests to shift the already existing kink to the nearest point in time, the corresponding variation of the trajectory is 
accepted with probability ~ 0(1). Thus, in average, it takes some l/(tAr) 2 attempts to create a new kink - anti-kink 
pair, and l/(tAr) attempts to move the kink to the nearest position in time. To conclude, the algorithm is becoming 
progressively ineffective when At is reduced. 

In this paper we prove that for any Hamiltonian with discrete Hilbert space it is possible to arrange an exact MC 
process, i.e., when both the description of the trajectories in time and their actions contain no approximations at all, 
and which is at least 1/(£At) 2 times faster than the conventional scheme in accumulating statistics. 

First, we observe that for discrete systems like (||) any trajectory may be easily described using continuous time 
variable: it is sufficient to define the state at t = and finite number of points in time {ti, T2, . . . , t„} where the state 
has changed from some \a>k-i) to |a/c). Within this description the action of the trajectory is known exactly, since 
one may formally think about taking the limit At — * while calculating S within the standard approach. When we 
turn to the problem of generating new kink-antikink pairs, we formally face an incontestable obstacle: the probability 
to accept the trajectory with one extra kink-antikink pair is formally zero, W new /W id ~ (tAr) 2 — > 0, so that the 
algorithm does not allow to change the number of kinks at the level of comparing two particular trajectories. This 
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obstacle, however, can be circumvented, if one knows an exact integrated statistics of trajectories with a given number 
of kinks. One may decide then about the number of kinks on the trajectory by relating these integrated probabilities. 
In what follows we demonstrate explicitly how to find this integrated statistics for arbitrary discrete system in any 
time window of width r c by relating it to the statistics of the biased two-level system. 

2. Let Hq and V be the diagonal and non-diagonal parts of the Hamiltonian, fl, in a chosen representation 
corresponding to the full set, {a}, of eigenstates of H , and H | a) — E a \ a). Then statistical operator can be 
ordinarily related to the Matsubara evolution operator, a, in the interaction representation, i.e., we write e~@ H = 
e~ 0H °u with 

<r=l- f drV(r) + ... + (-ir f dr m --- P dr x V{r m ) ■ ■ ■ V{n) + . . . , (3) 
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where V(t) — e TH °Ve~ TH ° . Without loss of generality and in accordance with typical forms of Hamiltonians of 
interest, V can be written as a sum of elementary terms, Q s , whose action on any function from the set {a} results 
in another function from this set: 

V = Y,Q>\ Qs\cc)=-q ia {s)\ 1 ) ( 7 = 7 ( S ,a)). (4) 

s 

Since V is Hermitian, for any s in the sum (0) there exists s' such that Q s i — Q\. Then we may rewrite Eq. (0) in 
components (below E ai = E a — E^): 
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Note that there is no additional summation over the indexes of the intermediate complete sets (labeled by Greek 
letters) since these are defined in a unique way by configurations of (si, s%, . . . , s rn ). 

We confine ourselves to the case of finite-range interaction, which is defined by the requirement that for each term 
si of elementary operators {Q s } there exists only a finite number of terms S2 for which the condition 

[Q Sl (n) ) Q S2 (r 2 )]=0 (6) 

is not met. In the case of finite-range interaction the structure of the series (g) is drastically simplified, the simplifi- 
cation being of crucial importance for practical realization of our algorithm. From (^|) it follows that up to irrelevant 
change in indexing energies and matrix elements, one may pay no attention to the chronological order of Q Sl (ti) and 
Qs 2 { T 2) in the evolution operator. This suggests representing of a general term of the series (^) in the following form. 
First, we introduce the notion of a "kink of kind s", which is characterized by a time r, a matrix element g Q7 (s), and 
a diagonal-energy difference E ai . The former two we will refer to as parameters of the kink. It is essential that (i) to 
obtain parameters of a kink one need not know explicitly the whole state | a), or | 7) - local information is enough; 
(ii) to set a particular structure of a term in Eq. (||), including the chronological order of all noncommuting operators, 
it is enough to specify for each kink its associated neighbors, i.e. nearest in time noncommuting kinks. 

Now we are ready to introduce a stochastic process directly evaluating Eq. (||). For simplicity, we assume that all 
q a /3(s) are positive real numbers. [In many particular alternative cases a straightforward generalization is possible, 
but usually at the expense of convergence.] Summations and integrations in Eq. (||) then can be regarded, up to 
a normalizing factor, as an averaging over a statistics of different configurations of kinks, each configuration being 
defined by a certain number of kinks of certain kinds, their associations and particular positions in imaginary time. 
Our aim is to organize a stochastic process of examining this statistics by generating different kink configurations in 
accordance with their weights. The process will consist of a number of independent sub-processes, responsible for 
modifications of particular types. The simplest sub-process is the random shift of the position of a kink in imaginary 
time without interchange with associated neighbors. The weight function for this sub-process is just exp(ri^ Q7 ), r 
is the position of the kink, E ai is the corresponding energy difference. It is clear that this process is far from being 
complete since it does not change neither the number of kinks nor their associations. 

The sub-process which is capable of doing so, and which will play the key part in the whole procedure, is the process 
of creation/annihilation of a kink - anti-kink pair. First we define the notion of the appropriate window (with respect 
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to kink - anti-kink pairs of the kind s - s'), as the interval of imaginary time of some chosen length To where (i) there 
is either one and only one kink - anti-kink pair of the kind s-s' (case A), or no such pairs at all (case -B), and (ii) 
there are no kinks associated with any kink of the pair. Now we notice that the total contribution of configurations of 
the class A (integrated over the positions of the pair inside the window and integrated and summed over all possible 
configurations of the rest of the kink pattern) differs from the total contribution of the configurations of the class B 
by the factor 



lg}{s) = J dr 2 J dT igai ( S ') e ^-g 7Q ( S )e T ' E -, (7) 
where r Q - is a randomly chosen left boundary of the time window. Elementary integration yields 
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This allows us to weigh these two classes. Having found an appropriate window we renew its contents by the following 
rule. First, we delete a pair s - s', if there is one. Then with the probability wa we create another kink - anti-kink 
pair s-s', and with the probability wb = 1 — wa leave it as it is. 

w A /w B = . (9) 

In the former case we also randomly choose the positions of the kink and anti-kink, in accordance with the distribution 
following from their exponential factors e' T2_Tl ' £ °T. 

It is worthnoting that the kink - anti-kink procedure actually relies on an important general fact, easily seen from 
Eq. (|]), that the local statistics of kink - anti-kink chains, corresponding to virtual transitions between two states, 
| a ) and | 7 ) , is exactly the same as for the biased two-level system with the bias energy equal to E aj and tunneling 
amplitude equal to q ai {s). The unnormalized probabilities of having exactly n kink-anti-kink pairs in the appropriate 
window are given then by 

r(2n)/^ _ (I ^( S ) I" T oT V* J n-l/2^) ~ ^+1/2^) ( 

a ^ {S) ~ n\ 2 6 (2x)«"V2 ' [W) 

where x = To-E a7 /2, and I u {x) is the modified Bessel function. Therefore, in principle a more general weighting 
procedure is possible, involving transformations of the kink - anti-kink chains of arbitrarily large length, which first 
deletes the whole chain of the kink - anti-kink pairs s-s' from the selected appropriate window and then creates 
another chain according to Eq. (|l(i|) . This, however, can hardly have practical importance because of small probability 
to find an appropriate window of large length. 

Generally speaking, the reason why the kink - anti-kink process in certain cases turns out to be incomplete is that 
there may exist more than one way to get from a state | a) to a state | 7). Consider, for example, a typical situation 
corresponding to the following two kink chains 

q 1 x{32) e T2E ^ q Xa (sx) e TlE ^ (case C) , (11) 

g 7 „(s 4 ) e TiE -<» q ua (s 3 ) e TaE - (case D) . (12) 

To obtain relative weights of classes C and D one could proceed in direct analogy with the kink - anti-kink processes, 
comparing the integral contributions. It is more reasonable, however, to take advantage of the fact that the number of 
kinks is the same in both cases. This allows us to take the differential limit. Namely, we compare only configurations 
when T\ and T3 {j2 and T4) vary in an infinitesimal interval dr around n = T3 = r a (T2 = T4 = T(,). As a result, we 
come to a Metropolis-type discrete sub-process of replacing the pair (si, S2) by the pair (S3, 84) and vice versa with 
probabilities wd and wc (respectively), satisfying 

w D q 1 u{s i )q va (s z ) 

= -, s -7—r exp (r b - T a )E X u } ■ (13) 

wc q-f\(S2)q\a(Sl) 

All the above-considered sub-processes deal with the "bulk" of the kink configuration. Finally, we introduce an 
"edge" sub-process, which, on one hand, changes one of the indexes of cr Q7 (for defmiteness, 7 into A), and, on another 
hand, changes the number of kinks by one. The appropriate window for the kink s is now defined as an interval [0, To], 
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where there is only one kink (case E) or zero (case F) kinks of the kind s, and no kinks associated with s . In close 
analogy with the kink - anti-kink sub-process, in the new configuration with the probability we there will be a kink, 
and with the probability wp — 1 — we it will be absent 



The position of the kink in the former case is chosen in accordance with the factor e T ^ A . 

To illustrate the above-described approach, let us formulate it for bosonic Hubbard models [and closely related 
to them (by Holstein-Primakoff transformation) spin lattices]. In the site representation, the diagonal part of the 
Hamiltonian, Hq, involves interparticle interactions and an external potential, if any. The nondiagonal part, V, 
corresponds to the intersite hopping, and terms Q s describe hopping from one particular site to another. The most 
transparent case is a hard-core one- dimensional system, especially if the quantities of interest are just potential and 
kinetic energies. As is easily seen, in this case the kink - anti-kink sub-process is enough to reproduce the statistics. To 
enhance the convergence, however, it is useful to utilize also the sub-process of the random shift of the kink position 
and the discrete sub-process of the kink - anti-kink pair reversal, i.e., Eq. ( |i"3[ ) with s 2 — s ii s 4 — s ii s 3 = s 'i- In 
d > 2 the discrete two-kink sub-process Eq. ( |l3| ) is essential to reproduce all possible trajectories of a particle. If one 
is interested in non-local (both in space and imaginary time) off-diagonal correlators or, say, the system behavior in 
a ring geometry with a global gauge phase, he should take advantage of the edge sub-process Eq. (M). 

It should be noted that the requirement that the interaction be of finite-range character is not of fundamental 
nature. In principle our process is applicable to interactions of arbitrarily long range, and, in particular, may work 
in momentum representation - if this is crucial for this or that reason. The only limitation which arises here is the 
small size of the appropriate window To. 

We note also, that our approach might prove to be the most effective and very accurate method for MC simulations 
of continuous systems if one first approximate them by lattice models. For one thing, in the case of harmonic oscillator 
the corrections for discreteness are known to be extremely small. 

3. All the above mentioned ideas were realized in a code allowing calculations of the one-particle properties for the 
Hamiltonian (0) . The results were compared with those obtained by the standard Monte Carlo world-line algorithm, 
based on finite At We found that, e.g., for the periodic potential t/j = (— l) l U and discrete harmonic oscillator 
Ui = Ui 2 the discrete-time code with At = 0.1/ 1 was more then two orders slower in accumulating statistics (or 
required > 10 2 longer computation time in order to achieve the same accuracy) than ours. This result is hardly 
surprising, because, as discussed above, the "bottleneck" in the standard procedure is connected with the unlikely 
process of the kink-anti-kink pair creation (with probability (tAr) 2 ~ 1CP 2 ). No matter how simple are our test 
models their description, nevertheless, contains all the essential features of simulation of much more complicated 
objects. 
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